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The photoabsorption spectra of calcium-doped argon clusters CaAr„ are investigated at thermal 
equilibrium using a variety of theoretical and numerical tools. The influence of temperature on the 
absorption spectra is estimated using the quantum superposition method for a variety of cluster 
sizes in the range 6 < n < 146. At the harmonic level of approximation, the absorption intensity is 
calculated through an extension of the Gaussian theory by Wadi and PoUak [J. Chem. Phys. 110, 
11890 (1999)]. This theory is tested on simple, few-atom systems in both the classical and quantum 
regimes for which highly accurate Monte Carlo data can be obtained. By incorporating quantum 
anharmonic corrections to the partition functions and respective weights of the isomers, we show 
that the superposition method can correctly describe the finite-temperature spectroscopic properties 
of CaAr„ systems. The use of the absorption spectrum as a possible probe of isomerization or phase 
changes in the argon cluster is discussed at the light of finite-size effects. 



I. INTRODUCTION 

Large van der Waals clusters offer a convenient medium 
to study solvation at the molecular level, and also to 
undertake chemical reactions with well defined thermo- 
dynamical conditions .ElcrB Reactivity of chromophore- 

\ doped rare-gas clusteijsJma been the subject of intense 
experimental effort .□eIeIoEIu Photodissociation of HBr on 

] argon-clusters has been investigated by Baumfalk et 
aLpffl including a theoretical study of the HBr-cluster 
interaction.cl Chemiluminescence in the bimolecular reac- 
tion Ba-I-N20 p-^BaO*-l-N2 was the subject of a work by 
Gaveau et al.p and photodetachment in IHI~Ar„ clus- 
ters was used to probe transition state spectroscopy of 
the I+HI reaction by Liu and coworkers.^ One may also 
cite the photoinduced reaction dynamics between Ca* 
and HBr towards CaBr*-|-H at the surface of a large ar- 
gon cluster carried out by Briant et alu Very recently, the 
reaction between N2O and Li2 mediated by Ar„ clusteLS 
has been investigated by Breckenridge and coworkers £j 

A practical understanding of the physical and chcmi- 
' cal properties of the deposited species relies on a good 
characterization of their spectroscopy. In turn, spec- 
troscopy gives great insight into the structure of the van 
der Waals complex, especi ally jp . concern with the lo- 
cation of the chromophoreJlillljIlJ This idea was used 
to propose a way of deteeting structural transforma- 
tions inside the clusterMla Not long after its was dis- 
\ covered in simulations,t£l the so-called dynamical co- 
existence phenomenon between solidlike and liquidlike 
phases received some experimental evidence by Hahn 
and Whetten who reported changes in the excitation 
spectrum of |J3jenzene-(argon)„ clusters with increasing 
temperaturcJlj Similar experimental results were re- 
ported for dichloroanthracene-Ar„ clusterS|-liy-| Even et 
aZ.tj More recently, Curotto and coworkerst3ll3'E£l theo- 
retically investigated the infrared spectroscopy in Ar„- 
HF clusters. They observed significant variations in the 
red shift associated with the two lowest vibrational states 



of the HF molecule, which they could pspsjipiate with the 
thcmiodynamical state of the cluster .tZHl§Ej Moseler et 
a/.EJ found that the photoabsorption spectra of Nag dis- 
plays some features characteristic of isomerization and 
structural fluctuations. Lastly, simulations on charged 
rare-gas clusters have investigated theJnfluence of tem- 
perature on the absorption spectrum. Ej'Ej So far, no di- 
rect evidence for isomerization in such rare-gas clusters 
has been found yet. 

While spectroscopy has many attractive features for 
experimental studies in large clusters, several difficul- 
ties pave the way for efficient theoretical approaches. 
Firstly, because of the numerous degrees of freedom, a 
fully quantum calculation of the absorption spectrum by 
exact methods (wavepacket or diagonalization in a ba- 
sis) is not feasible beyond a few atoms and one is usually 
lead to simplify the problem. A common approxima- 
tion is the Condon hypothesis, which assumes vertical 
transitions. Secondly, the equilibrium properties at fi- 
nite temperature require some sampling of the ground 
state potential energy surface (PES). This sampling can 
be hard to get because of the topography of the multidi- 
mensional PES (the so-called energy landscape). When 
the energy landscape is complicated or has multiple fun- 
nels, the relaxation to equilibrima may be particularly 
slow with a glassy-like behavior .E^ Broken ergodicity is 
then likely to occur in conventional molecular dynam- 
ics (MD) or Monte Carlo (MC) simulations, and also in 
quantum simulations. Fortunately some remedies to this 
problem have recently been piggested, including multt. 
canonical ensemble samplingHJ or parallel tempering.Lj 
Thirdly, most experiments involving chromophore-doped 
rare-gas clusters take place at low temperatures, at which 
the classical approximation may be questioned, especially 
for the lighter elements neon and helium. Path-integral 
Monte Carlo or centroid MD should then be used instead 
of classical techniques, however their use with ergodicity 
restoring strategies remains numerically demanding for 
systems containing tens or hundreds degrees of freedom. 
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We have recently proposed an alternative way of 
calculating thermodynamical observables at canonical 
equilibrium, using tlie superposition approximation of 
quantum oscillators .Ej This method, originally devel- 
oped for thermodynamical observables only, was im- 
proved ph cnQrH cnologically by including anharmonic 
corrections J23EZI and extended to generally treat other 
temperature-dependent observables .Ej Its main interest 
is to circumvent the troubles of broken ergodicity by con- 
sidering separately the contribution of each isomer to the 
partition function. The method also naturally describes 
quantum delocalization in a straightforward way. 

One goal of this paper is to extend the quantum su- 
perposition method to the realm of calculating photoab- 
sorption spectra. In addition to the ingredients presently 
available, one needs to calculate the photoabsorption 
spectrua of each stable isonaer. For this, we use a recent 
theory by Wadi and PoUakcJ that we have extended to 
the general case of arbitrary potential surfaces. As will be 
seen below, this Gaussian theory of absorption succeeds 
in predicting absorption line shapes, widths and shifts for 
harmonic polyatomic molecules in good agreement with 
highly accurate Monte Carlo data. 

With the quantum superposition method, we are able 
to carry out systematic studies for a given set of clus- 
ter sizes. We can then address our second goal, and try 
to relate the features of the absorption spectrum to the 
knowledge of isomerizations and possible phase changes 
inside the cluster. More specifically, by looking at par- 
ticular sizes we hope to gain insight in the relationship 
between structure, temperature, quantum delocalization, 
and the shape and position of the spectral lines. 

This work follows previous efforts aimed at inter- 
preting ejiperimental absorption spectra on large CaAr„ 
clusters.t2l The spectroscopy of the CaAr diatomic 
molecule has been investigated both experimentally and 
with ab initio calculations in Ref. and a Diatomic- 
In-Molecules (DIM) Hamiltonian was constructed to 
model the. 4s ground- and 4s4p excited-states of CaAr„ 
clusters.Ej In the present context of finite-temperature 
spectroscopy, the heavy numerical cost of statistical sam- 
pling is balanced with this relatively cheap, albeit ac- 
curate Hamiltonian. The paper is organized as follows. 
In the next section, we briefly recall the quantum su- 
perposition method in the harmonic approximation. We 
also include the derivation of approximate perturbative 
anharmonic quantum corrections to the partition func- 
tions. The Gaussian theory of absorption for hamionic 
molecules, initially proposed by Wadi and Pollak,cj will 
be generalized to include possible Duschinskii rotations 
in Sec. [II. These methods will be tested on simple clus- 
ters, in both the classical and quantum regimes, by com- 
paring their results to high statistics Monte Carlo data. 
Their main application to the variations of the photoab- 
sorption spectrum with temperature in some size-specific 
CaAr„ clusters will be presented in Sec. iV. The results 



the paper in Sec. ^ 



II. QUANTUM SUPERPOSITION METHOD 

In the superposition approach EElii also known as 
the multiple normal modes modelj^ the total partition 
function Z of the iV— atom system in equilibrium at tem- 
perature T = l/ksP is replaced by a sum over stable 
minima {a}: 



(1) 



Each minimum a, or inherent structure, has a contribu- 
tion Za to Z, further weighted by the degeneracy factor 
na, which accounts for permutational isomers and global 
symmetry. In CaAr„ clusters {N = 1), is given by 
2n\/ha, where ha is the order of the point group of struc- 
ture a. The problem has now been reduced to estimating 
the individual functions Za- At the simplest level of ap- 
proximation, the harmonic vibrational frequencies {^0,^} 
provide a zero-th order expression for Za'- 
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with h the reduced Planck's constant and Ea the po- 
tential energy minimum of isomer a. In the classical 
regime h ^ 0, the decomposition (|l|) is formally exact, 
since each point of the configuration space (except saddle 
points) belong to one and only one basin a. Lo-the quan- 
tum regime (finite K), tunnelling is neglected.E3 In princi- 
ple, the knowledge of all minima of the energy landscape 
allows us to calculate any thermodynamical observable, 
including internal energy and heat capacity, from Eq. (Q). 
However, for all but the smallest atomic systems, exhaus- 
tive enumeration of these isomers is not possible because 
the number of local minima is thought to grow i^spo- 
nentially with the number of degrees of freedom.tJ To 
correct for the uncompleteness of the available set of iso- 
mers, a reweighting technique can be usedBjS The idea 
is to perform a simulation at a reference temperature Tq, 
and to compute the probability Pa of observing isomer 
a- The partition function Z is approximated as a partial 
sum over weighted isomers: 
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with the partial character of the sum indicated by the 
prime symbol. If we assume that the weights {ga} are 
independent of temperature, the probability Pq(T'o) of 
finding isomer a at Tq is proportional to ganaZa{Po)- 
Hence, Z can be expressed as 



will be discussed in view of the isomerization and phase 
changes phenomena. Some concluding remarks will close 
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A classical simulation can be used, to sample the min- 
ima and their respective weights. l3 In the above equa- 
tion, the denominators Za{Po) should then be replaced 
by their classical analogue. Even though the resulting 
partition function (Q) should not depend on the classical 
or quantum nature of the simulation, some differences 
can be seen when using semiclassical simulations includ- 
ing quantum corrections^ This reweighting procedure 
has been shown to improve significantly the predictions 
of the superposition method. In particular, Wales was 
able to obtain a 'S-bend' in the kinetic temperature of 
Ar55,t3 as seen previously in simulations. t3 

The other important limitation of the method is the 
harmonic approximation used to calculate Za{P). Sev- 
eral phenomenological anharmonic corrections hajvepbeen 
proposed in the litterature for classical systemsOO In 
the case of quantum systems, the optimal normal modes 
method of Cao and VothEa can be used to include anhar- 
monicity to some extent, but the self-consistency required 
makes it computationally heavy for large sets of isomers. 
Building upon a recent effort to include systematic cor- 
rections from the enpgy surface topography via classical 
perturbation theoryJ23 we can attempt to include anhar- 
monicities in the quantum partition functions Za in a 
more rigorous fashion. Since the classical corrections are 
also numerically demanding, we restrict ourselves to low- 
order corrective terms. In order to make the reading eas- 
ier, the details of the perturbation expansion are given in 
the Appendix. The superposition method not only pro- 
vides the thermodynamical equilibrium properties, and 
other quantities can be computed-provided they can be 
calculated for each minimum aEB At thermal equilib- 
rium, the average value (A) of the global observable A is 
the weighted sum over minima: 
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III. GAUSSIAN THEORY OF ABSORPTION 

We assume that the cluster is in canonical equilibrium 
at temperature T = l/fc^/?, and that it undergoes an ex- 
citation by an infinitely narrow pulse at frequency lu from 
its ground state potential surface V^(R). The normalized 
excitation spectrum isH 



X(c^,/3) = 



E„ 



with \gi) and \ej) the eigenstates with eigenenergies , 
Ef,- of the ground- and excited-state Hamiltonians Hg 
and i7e, respectively. The transition dipole moment 
Mij = l(6jlAi|<?i)P will be taken as constant in the usual 
Condon approximation of vertical excitations. is the 
energy difference between the bottom of the ground- and 
excited-state surfaces. The absorption intensity can be 



expressed as the Foucipr transform of the imaginary-time 
correlation functioni 
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x(r,/5)e-^^(^^-''-"'dr, (7) 



with 
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and where Z{(3) — Tre^''^^a-,is the partition function. 
Following Wadi and Pollak,L£l we now perform a short- 
time expansion in the correlation function. This expan- 
sion should be valid at moderate temperatures and for a 
large number of degrees of freedom. Indeed, in such cases 
the combination of the numerous vibrational modes will 
significantly dephase rapidly, and the correlation func- 
tion will be nonzero only for short times. Using this ap- 
proximatio]x,.jthe absorption intensity can be calculated 
analyticalljo as a Gaussian function: 
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with ujQ — Eq/H, Eq being the vertical transition energy. 
The line shift Aoj and the width w are given as a func- 
tion of the two first moments of tho-ahiftcd Hamiltonians 
difference AH = He- Hg- Eq byfea 



hAuj = (AH), 



and 



hio = (AH' 
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We have noted here (A) the thermal average of observable 
A over the ground state potential surface Hg. We now 
write the Hamiltonians Hg and He in harmonic (normal 
modes) coordinates Q: 



H, 



He = 



1 



(12) 



Pi 



GtQ + iqtWeQ + So. (13) 



In these equations, is the mass- weighted Hessian ma- 
trix of the ground state at the lowest energy minimum 
Qo, which is assumed to be located at the origin of the 
reference frame. Wg can be written as = diag ('j-'g. ), 
as it is diagonal by construction. We = [w]^) is the sec- 
ond derivatives matrix of the excited state surface at Qq. 
Lastly, G = 9_ffe(Qo)/9Q is the excited state gradient at 
Qo. The matrix We is not assumed to be diagonal, and 
may include the so-called Duschinskii rotations.cil It is 
then a simple matter of algebra to obtain the expressions 
for Ao; and Q. Using the notation pi = LUg.ta'Dh(3hujg-/2, 
we find 



Auj 
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Apart from the different notations, the above expressions 
differ from the results by Wadi and PoUak by the second 
term in the right-hand side of Eq. (|l^), which was ne- 
glected by these authors Although Eqn. (|lj) and (|l5| ) 
are only valid for polyatomic, harmonic systems, anhar- 
monic corrections can be carried out in the very same 
way as the perturbative expansion mentioned in the pre- 
vious section and detailed in the Appendix. However, 
the present expressions are aheady the outcome of an 
approximation in the correlation function, hence we did 
not attempt to include such corrections here. In addi- 
tion, their estimation would require the computation of 
higher order derivatives of the excited states energy sur- 
faces, which will be a heavy numerical task, even for for 
DIM potentials. 

We have tested the Gaussian theory on small CaAr„ 
clusters. The ground- and excited-state Hamiltonians 
have been fully described previously.tSi They are-based 
on the Diatomic-In-Molecules (DIM) formalism,c3 and 
were parameterized on ab initio calculations .e3 To give a 
better account of experimental properties of the CaAr di- 
atomics, we have slightly modifietLthe parameters with 
respect to our previous works.ESEa The Ca-Ar ground 
state pairwise potential has been scaled by 0.85, and the 
excited state surfaces have been shifted by -1-305 cm^^ 
to fit the 4s4p experimental transition in calcium. As a 
consequence, some global minima may have changed, as 
in CaAray. 

In the present work, we focus on the 4s^ AsAp 
transitionii^ which have received some experimental 
attentionJlj The DIM Hamiltonian provides the energy 
and numerical derivatives, as well as the electronic tran- 
sition dipole moment, for any atomic configuration R. It 
can thus be used to simulate absolute absorption spectra. 
The ground state vibrational eigenmodes and eigenfre- 
quencies {lOq^ } were obtained from the diagonalization of 
the Hessian matrix whose elements were computed from 
their analytical expressions. In both the classical and 
quantum cases, the absolute spectrum is calculated as 
the sum over the three possible excited states. For the 
Gaussian theory, the three gaussian curves are weighted 
by the square of the transition dipole moment at the 
equilibrium ground state geometry. 

First, we have calculated the classical spectrum from 
the Gaussian theory in the classical limit ?i — > 0, and from 
classical Metropolis Monte Carlo simulations at low tem- 
perature, using 10^ cycles following 2 x 10^ equilibration 
cycles. The spectra of CaAris and CaAr54 at T = 5 K are 
displayed in Fig. |l|. For comparison, we have also plotted 
the predictions of the Gaussian theory in the classical 
limit ?i — > 0. In this limit, the line shift and line width 
tend to their classical values (7iAa;)c and {hu))1 given re- 
spectively by 

(hAu), = , (16) 




The reliability of our comparison was checked by ensuring 
that only a single isomer was visited during the Metropo- 
lis MC simulation. As can be seen in Fig. |l|, the accuracy 
of the Gaussian theory is very good at such low tempera- 
tures. The global minimum geometry of GaAri3 has two 
spectral lines located at 23658 (2 degenerate lines) and 
24034 cm~^, respectively. Due to its finite temperature, 
the cluster loses its instantaneous symmetry and the de- 
generacy is lifted. This explains why three distinct lines 
are seen in Fig. |^(a). In GaAr54 the Gaussian theory of 
absorption remains in good agreement with Monte Garlo 
data. This cluster has again 2 degenerate sharp lines at 
23715 cm-\ and one broad hue at 24080 cm"!. This 
time, the degeneracy is no longer lifted by temperature, 
and the two lines still merge into a single peak. In the 
classical theory the absorption intensity can also be ob- 
tained analytically by integrating Eq. (^ under the Con- 
don assumption without performing the short-time ex- 
pansion. Therefore, in the classical case, the differences 
between the calculated and simulated spectra in Fig. |l| 
can be attributed to anharmonicities only. 

We now turn to a more complicated situation, namely 
CaAr37 in the quantum regime, and at moderate temper- 
ature T = 25 K. We will show below that this tempera- 
ture is located just above the melting point of this cluster. 
Many isomers are present, and this will provide a test 
of the superposition approximation. In the simulation, 
quantum effects have been included in a quasiclassical, 
corrective way byi-adding the corrective Feynman-Hibbs 
effective potential to the ground state energy. This po- 
tential was shown to give good results on argon clusters 
with respect to the more accurate Fourier-path-integral 
Monte Carlo simulation of Neirotti and coworkers,a es- 
pecially when compared with the alternative Wigner- 
Kirkwood corrections .l3 Since it is quasiclassical in na- 
ture, it cannot be used in the low temperature regime 
T ^ 5 K where only few isomers are present, and path- 
integral simulations would have to be carried out in- 
stead. Because the superposition method does not as- 
sume any barrier nor pathway between the isomers, it 
yields supposedly ergodic data within the set of available 
minima. Therefore a reliable comparison with simula- 
tion results must also correct for possible broken ergod- 
icity, and for this we have used parallel tempering Monte 
Carlo.El Also, since the effective potential explicitely de- 
pends on temperature, the probability of accepting an 
exchange between, adjacent trajectories should be modi- 
fied accordingly.^ 

Concerning the superposition approximation, and to 
quantify the importance of perturbative corrections 
in the partition function, we have separately consid- 
ered zero-th order harmonic as well as second-order. 
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FIG. 1: Classical absorption spectra of the 4s^ 4s4p transition in (a) CaAna and (b) CaAr54 from Monte Carlo simulations 
(dashed lines) and the Gaussian theory (solid lines) at T = 5 K. 
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FIG. 2: Absorption spectra of CaArsr at T = 25 K from 
quasiclassical Monte Carlo simulations (dashed line) and the 
quantum superposition method in the harmonic approxima- 
tion (long dashed line) or with anharmonic perturbative cor- 
rections (solid line). 



anharmonic-corrected partition functions with the same 
database of isomers. The absorption spectra obtained 
with the Gaussian theory in the superposition approach 
and with quasiclassical MC simulations are represented 
in Fig. H. The lowest energy isomer of this cluster, includ- 
ing the zero-point energy contribution, has three distinct 
spectral lines for the 4s^ AsAp transition, namely at 
23840, 23970, and 24260 cm^^. As can be seen from 
Fig. g, none of these lines remain visible as peaks. We 
observe a very good agreement between the present the- 
ory and on-the-fly simulation, indicating that the relative 
weights of the isomers are correct. However, at this point. 



we cannot tell whether the visible shifts and broadenings 
are an effect of temperature only, or of the different spec- 
troscopic signatures of the other isomers present. Some 
shifts in the peaks remain slightly approximate in the 
superposition approach, and the second-order corrected 
weights only marginally improve the results. At 25 K, 
quantum effects turn out to be rather small on the line 
width and shape, but they may be more important on 
the weights of the isomers, as will be seen below. There- 
fore, from Fig ^ we see that anharmonic corrections to 
the isomers weights have a relatively weak influence for 
quasiclassical clusters. More information on the impor- 
tance of these corrections will be gained by looking at 
the thermodynamical curves in the next section. 

The Gaussian theory presented in this section can be 
extended to purely quantum systems at T = 0, for which 
delocalization is only due to the quantum nature of the 
vibrational modes, but not to an actual atomic displace- 
ment induced by a finite energy deposit. In this pur- 
pose, the Boltzmann weights in Eqn. (^ and (|^) are re- 
placed by the probability density of finding the system in 
its ground state at configuration R as the square of the 
wavefunction. In the harmonic approximation, this prob- 
ability density is a multidimensional Gaussian function, 
and the integrated absorption spectrum is found to be 
Gaussian with line shape and width given by the /3 — > oo 
limits of Eqn. (pll) and (p^): 
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(a) Quantum CaAr^ T=0 K 



(b) Quantum CaAr,3T=0 K 
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FIG. 3: Absorption spectra of (a) CaAr2 and (b) CaAna at T = K from diffusion Monte Carlo simulations (solid lines) and 
the Gaussian theory (dashed lines). For CaAr2, we have also represented the results of wavepacket relaxation calculations (long 
dashed line). 



The behavior of the Gaussian theory in the T = quan- 
tum case has been tested on CaAr2 and CaAri3 using 
simple diffusion Monte Carlo (DMC) in internal coor- 
dinates (without global rotation and translation), and 
without importance sampling. The DMC calculations 
were carried out with 1000 replicas, and the wave func- 
tions were averaged over 10^ steps after 10^ initial equili- 
bration steps. In addition, for CaAr2, the vibrational 
ground state wavefunction was also computed numer- 
ically on grids, using quantum wavepacket relaxation 
(WPR) method (imaginary time propagation) .cZl The 
propagation was performed for zero total angular mo- 
mentum, representing the wavepacket by an harmonic os- 
cillator Discrete Variable Representation (DVR) for each 
of the two Ca~Ar bond coordinates, and by a Legen- 
dre DVR for the bending angle ArCaAr. The primitive 
basis for representing one bond coordinate includes 10 
harmonic oscillator eigenfunctions. Associated Legendre 
functions with j < 60 and m — were used for repre- 
senting the bending mode. In both the DMC and WPR 
methods, which are both supposed to give essentially ex- 
act results, the absorption spectrum was calculated as- 
suming vertical transitions from the quantum sample of 
geometries. 

In Fig. ^ the absorption spectra calculated with these 
various methods are compared. For the small, triatomic 
cluster, the three methods give very similar peaks. In 
its equilibrium geometry, CaAr2 has three nondegener- 
ate spectroscopic lines at 23621, 23670, and 23918 cm~^, 
respectively. Quantum delocalization of the vibrational 
modes leads to a rather small shift of these lines, and to a 
significant broadening. For CaAri3, the simulated peaks 
are also in good agreement with the Gaussian prediction, 
but we notice that some asymmetry is visible on the "ex- 
act" spectrum from diffusion Monte Carlo. These effects 



are clearly beyond a Gaussian description, and are due to 
anharmonic effects. These effects are relatively weak for 
the rather heavy calcium-doped argon clusters but may 
become more important for neon or helium clusters. 

The above results suggest that the Gaussian thfe 
ory of photoabsorption developed by Wadi and Pollakn^ 
upon the formal results of Yan and Mukamelcj can give 
quite accurate results in both the classical and quantum 
regimes at low and moderate temperatures. 



IV. SPECTROSCOPIC SIGNATURE OF PHASE 
CHANGES 

The Gaussian theory described in the previous Section 
will now be used with the superposition approach to in- 
vestigate the influence of temperature on the absorption 
spectra in context with isomerization phenomena. 

For the two smallest clusters chosen, CaArg and 
CaArio, we performed extensive sampling of the ground 
state PES, and we believe that nearly all stable isomers 
have been identified. To apply the superposition method 
to larger CaAr„ clusters, the reweighting procedure from 
simulations must be used. Quasiclassical parallel temper- 
ing Monte Carlo simulations of 10^ + 2 x 10^ equilibra- 
tion cycles were used, for a set of temperatures equally 
spaced by 2.5 K in the range 7.5 K< T < 60 K for 
n < 54, and 15 K< T < 75 K for CaArug. The iso- 
mers were found by quenching the 50 K trajectory for 
all sizes except CaAr37, for which the 25 K trajectory 
was used, and CaAri46, quenched at 60 K. Local mini- 
mizations were performed every 500 cycles. To prevent 
evaporation, a rimd spherical container was used in the 
MC simulation.E3c3 The radii of the containers and the 
number of different isomers obtained from the quenches 
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FIG. 4: Global minima of CaAr„ found by the basin-hopping algorithm after including the harmonic zero-point energy, (a) 
n = 6; (b) n = 10; (c) n = 13; (d) n — 30; (e) n — 37; (f) second lowest energy minimum of CaArar; (g) n — 54; and (h) 
n = 146. 



are given in Table | for the set of sizes chosen. We have 
also indicated the point group and energy of the global 
minima after including the zero-point harmonic contri- 
bution. These lowest-energy structures were generated 
from the basin-hopping or Monte Carlo + minimization 
algorithm.EJ They are shown in Fig. ^. In most cases, 
the global minimum of CaAr„ is very similar to the one 
of A r„4.i, with one argon atom substituted by calcium. 
Due to the lower binding energy andJarger equilibrium 
distance of CaAr with respect to Ar2j^ calcium is always 
found on the surface of the remaining argon cluster. The 
nearly identical masses of calcium and argon atoms also 
induces relatively small changes in the vibrational prop- 
erties, and the similarities between CaAr„ and Ar„+i 
remain when including the zero-point corrections. How- 
ever, in the case of CaAr37, we find the somewhat sur- 
prising result that the global minimum geometry is deca- 
hedral, and remains so after including the harmonic zero 
point energy contribution. Slightly higher in energy, sev- 
eral icosahedral structures are found, either in the form 
of uncomplete multilayer icosahedra (Mackay type) or as 
polyicosahedra (anti-Mackay type) . As in Arsg , quantum 
delocalization favors anti-Mackay structures over Mackay 



shapes]^ yet a Mackay structure is found as the sec- 
ond isomer in both the classical and quantum regimes. 
This isomer is represented in Fig. 0(f). Additionally, the 
distorted truncated octahedron geometry lies above, but 
very close, to these isomers. As a consequence, the en- 
ergy landscape of this particular cluster reveals interest- 
ing features, which will be investigated more deeply in the 
following article of this series .E2I This result is also consis- 
tent with the richer thermodynamical behavior generally 
observed in mixed van. dfir .Waals clusters with respect to 
homogeneous clusters. eHeHI 

Therefore, even though CaAr„ and Ar„_(_i clusters 
share many structural features, they also exhibit signif- 
icant and unexpected differences, which could affect the 
thermodynamics and dynamics in large extents. 



A. Isomerization in CaAre and CaArio 

The lowest energy minimum of CaArg is a distorted 
pentagonal bipyramid. The heat capacity of this clus- 
ter, computed from the superposition approach with or 
without anharmonic perturbative corrections, is repre- 
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TABLE I: Properties of the CaAr„ clusters simulated in this work. The energy and point group refer to the lowest energy 
minima, including the harmonic zero-point contribution. 
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sented in Fig. ^(a). For comparison, we have also plotted 
on the same curves the direct results of parallel tem- 
pering Mpjjite Carlo with the Feynman-Hibbs effective 
potential. EJ As is well known for this potential,c3 the 
canonical heat capacity diverges upwards at low tem- 
peratures. Due to its small size, this cluster shows a 
nearly featureless caloric curve, and a flat heat capacity 
(except the low temperature increase) in Fig. |(a). The 
small hump at T ~ 20-25 K is indicative of isomeriza- 
tions, and is followed by another bump near 45 K in the 
MC simulation. This latter broad peak is in fact spu- 
rious, as it comes from the dis^ewition of the weakly 
bound Ca atom. Previous worksSo have identified this 
peak as the signature of the liquid-vapor transition, from 
its dependence on the container radius. Here it would 
have been difficult to employ a radius smaller than 6A, 
because the motion of argon atoms would have been 
hindered. The 25 K hump is barely visible, because it 
is hidden and smoothed out by quantum effects. The 
agreement between the superposition approximation and 
the quasiclassical MC simulations is very good except for 
the extra liquid-vapor bump. This feature is obviously 
missing from the assumptions made in the superposition 
method. The quantitative effect of perturbative correc- 
tions is small for this cluster. 

CaArio has the ground state of an uncomplete icosa- 
hedron. Some isomers involving different locations of the 
calcium atom lie not far above the global minimum, re- 
sulting in a shoulder in the heat capacity near 10-15 K 
in Fig. |^(b). The main peak near 30 K is due to a more 
complete rearrangement of argon atoms, and is repro- 
duced quantitatively by the superposition approach after 
including perturbative corrections in the partition func- 
tion. 

The absorption spectra of these two clusters are rep- 
resented in Fig. 1(c) and 1(d) as a function of temper- 
ature. In each of these figures, the classical spectra ob- 
tained from the limit 7i — > has been superimposed to 
the quantum spectra. For the two sizes, the patterns are 
rather similar: the spectrum varies from a clear three- 
peak structure at low temperature towards a mostly sin- 
gle peak shape at high temperatures. The transition oc- 
curs roughly at 20 K in CaArg, and near 30 K in CaArio . 



In the case of CaArg the high temperature main peak 
appears made of two smaller peaks, in agreement with 
simulations .E3 From a general point of view, quantum ef- 
fects mainly consist in a broadening of the peaks. As 
expected, they get smaller and smaller as temperature 
increases, and becomes barely visible above 25 K. 

Therefore, in these two small clusters, absorption spec- 
troscopy seems to be an effective mean of detecting iso- 
merization events, as the range of variation in the spec- 
trum approximately matches the width of the heat capac- 
ity changes. As cluster size grows, one can thus expect 
to see sharper transitions, and to possibly use the cal- 
cium chromophore as a probe of the argon cluster melting 
point. 



B. Phase changes in CaAr„, 13 < n < 54 

We have chosen a set of sizes in the range 13 < n < 54, 
spanning between 1 and 2 icosahedral layers, and which 
are likely to exhibit significant finite size effects. We 
first discuss the thermodynamical behavior of these clus- 
ters. The heat capacities obtained using the superpo- 
sition method, with harmonic or second-order corrected 
anharmonic partition functions, are represented in Fig. ^ 
for the clusters corresponding to n = 13, 30, 37, and 
54. In the case of the larger cluster CaAr54, and also for 
CaAri46, it was not possible to include all the required 
anharmonic corrections due to the too large computa- 
tional requirements. Since they are not the main focus 
of the present paper, these thermodynamical data will be 
only briefly discussed in order to guide our interpretation 
of the following absorption spectra. 

In general, we find again a rather good agreement be- 
tween the quantum superposition results and the Monte 
Carlo simulations for all four sizes. As expected, anhar- 
monic corrections do have a positive quantitative effect 
by bringing the superposition data closer to the simu- 
lation results. However, no new qualitative feature is 
visible when they are taken into account. All heat ca- 
pacities exhibit a maiur^ak indicative of the solidlike- 
liquidlike phase change.La This peak is located near 40- 
45 K for CaAri3, CaArao and CaAr54, and near 20 K for 
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FIG. 5: Heat capacities computed from quasiclassical Monte Carlo simulations (full circles) and from the superposition method 
with harmonic (dashed lines) or anharmonic (solid lines) partition functions for (a) CaAre and (b) CaArio. Also represented 
are the finite-temperature absorption spectra calculated from the superposition method and the Gaussian theory of absorption 
in the quantum (solid lines) and classical (dashed lines) regimes for (c) CaArs and (d) CaArio. 



CaAr37. This lower melting point is similar to what is 
observed in the 38-atom ^ax^n cluster described by the 
Lennard-Jones potential.EjEll In addition to this main 
peak, extra humps or even peaks can be observed at 
lower temperatures in CaArso and CaAray. These fea- 
tures are indicative of preliminary isomerizatioi* result- 
ing from solid-solid transitions in these clustersEa Again, 
they bear similarities with the known finite temnerature 
behavior of LJ clusters with the same size.ES&OOCJ In 
CaAri3, the quantum character of the specific heat hides 
again a shoulder in the classical curve, as is the case in 
CaArg. This cluster happens to have a much larger num- 
ber of stable isomers than Ari4. In particular the capped 
icosahedron can be found with the calcium atom either 
in adatom location (lowest energy structure) or inside 
the icosahedral shell, resulting in a set of 4 different iso- 
mers slightly higher in energy. Structures not based on 



the capped icosahedron, or with the calcium atom in the 
center of the icosahedron, are less stable than these five 
geometries. Therefore the energy spectrum of local min- 
ima can be put into three sets, one being the global min- 
imum only, one being the other capped icosahedral iso- 
mers with calcium in the icosahedral shell, and all other 
isomers not build on the capped icosahedron motif. Be- 
cause of these three sets, the melting process involves es- 
sentially two steps. In Fig. ^(a) we plot the probabilities 
of finding the CaAria cluster in its ground state, or in one 
of the 4 subsequent isomers, or in the remaining, liquid 
configurations. These probabilities are computed using 
the anharmonic corrected quantum partition functions. 
They display a clear two-step isomerization pattern. In 
the first step, at about 30 K, the calcium atom is able to 
leave its capping site in favor of one of the twelve icosahe- 
dral sites. Then, above 40 K, more disordered structures 
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FIG. 6: Heat capacities of CaAr™ clusters, (a) n — 13; (b) n — 30; (c) n = 37; and (d) n — 54. The notations are the same as 
in Fig. |. 



become thermally available the melting peak takes place. 

In CaAr3o the preliminary solid-solid transition near 
2.5 K involves polyicpsaJiedral and multilayer icosahedral 
isomers, as in Ar3i.c3E3 Even though CaAr37 does not 
have the same global minimum as Arag, a similar ther- 
modynamic behavior occurs. We have calculated from 
the harmonic superposition approximation the probabil- 
ities of finding the cluster either in its decahedral ground 
state isomer, in its (Mackay-type) icosahedral isomers, or 
in any of the remaining isomers. The variation of these 
probabilities with increasing temperature are displayed 
in Fig. |^(b) . They show that the transition at about 5 K 
is due to a global decahedral-icosahedral rearrangement, 
and bears some similarity with the (C6o)i4 cluster .EiJ 
The anti-Mackay and truncated octahedral structures 
are only marginally populated at thermal equilibrium, 
and the main heat capacity peak is associated with the 
emergence of a large number of disordered configurations. 
Hence it can be considered as the signature of melting. 



For comparison, we have also plotted in Fig. |^(b) the cor- 
responding curves in the classical limit. While the global 
behavior remains essentially unchanged, the preliminary 
transition occurs at a higher temperature, near 10 K. 

The photoabsorption spectra of these clusters are rep- 
resented in Fig. ^for increasing temperatures in the range 
< T < 60 K for all sizes except CaAray, for which 
< T < 40 K. The classical spectra have also been 
plotted as dashed lines for T > 0. Except for CaAr37, 
the classical spectrum looks very much like its quantum 
counterpart with slightly shifted but noticeably narrower 
excitation peaks. The four clusters display different ab- 
sorption spectra, characteristic of finite-size effects. In 
CaAri3, the spectrum gradually evolves from a two-peak 
structure (at about 23700 cm~^ and 24100 cm^^) to 
another two-peak structure (at about 23650 and 23800 
cm"-'^). The peak ending at 23800 cm^^ appears first at 
30 K, and originates from isomers where the Ca atom 
lies in the main icosahedral shell. No real signature of 
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FIG. 7: (a) Probability to find the CaAria cluster with the calcium atom in a capping position (solid line), inside the icosahedral 
shell (dashed line), or in the liquidlike state (long dashed line), (b) Probability to find the CaArsr cluster in its decahedral 
ground state (black solid line), in a Mackay-type icosahedral isomer (black dashed line), or in any other isomer (black long 
dashed line). The corresponding classical curves are also plotted as grey lines. 



melting at 40 K is seen in Fig. ^(a), however the present 
spectra can still be interpreted in terms of isomerizations 
within the cluster. 

In the case of CaArao, the low temperature isomer- 
ization below 5 K has a strong effect on the absorption 
peaks, which are red shifted by about 200 wavenumbers. 
The spectra of this cluster also evolves very continuously, 
and the two main peaks at 23550 and 24000 cm~^ un- 
dergo a slight red shift and thermal broadening. The 
weak melting phase change is reflected on these spectra, 
and the most prominent feature remains the sharp iso- 
merization below 5 K. 

Near T = 10 K, the multiple peak structure of the 
spectrum in CaArsy also reveals that several isomers are 
coexisting with distinct spectroscopic signatures. The 
solid-solid transition seen in the heat capacity of CaArav 
at low temperature is reflected on the absorption spectra 
in Fig. ||(c). The three-peak T = K spectrum is pe- 
cuHar to the decahedral ground state. At T ~ 10 K, a 
significant blue shift occurs in the middle peak, provid- 
ing a signature of the icosahedral isomers. The quanti- 
tative difference in the solid-solid transition temperature 
induced by quantum effect is also seen in the absorption 
spectra between 5 and 15 K. Near 15 K, the spectrum 
can be clearly separated into more than three peaks, and 
suggests that distinct icosahedral isomers are present. At 
around 20 K, the various peaks merge into a double-peak 
structure, where the peak at 23800 cm~^ dominates the 
one at 24750 cm~^. Here we find a clear thermodynam- 
ical influence on the spectroscopic properties. Hence, in 
this cluster, the absorption spectrum contains notable 
informations about isomerization and phase changes, in- 
cluding some rather fine details on the characteristic tem- 
peratures. 



The equilibrium geometry of CaAr54 is a double layer 
Mackay icosahcdron where the Ca atom is located on 
a vertex site, thus keeping a relatively high symmetry 
(C5„). As in Ar55, the haal, capacity near peak 40 K 
indicates volume meltingJ^Cj while preliminary surface 
melting has a very small thermal signature. In Fig. ^(d) 
we observe that the two absorption peaks undergo a sig- 
nificant change near 30 K, and this corresponds to co- 
existing isomers with the Ca atom located on different 
substitutional sites of the external icosahedral shell. In 
this case the photoexcitation spectrum is sensitive to sur- 
face melting. At higher temperatures T > 40 K, many 
isomers are found in which the metal atom is more fully 
solvated by argon. The spectra show another qualitative 
change, as the peak near 23600 cm^^ disappears and is 
replaced by a broad bump near 24300 cm^^. 



C. CaAri46 

It is interesting to compare the situation in CaAr54 
with the one in the larger Mackay icosahcdron, CaAr 145. 
As in CaAr54, the ground state equilibrium geometry has 
C^y symmetry, with calcium located in a vertex site of the 
external icosahedral layer. This cluster has a huge num- 
ber of isomers, and many high lying minima are missed 
during the quenching process. In the vast majority of 
structures in our sample the calcium atom remains in 
the outer argon layers. This is not surprising, since, and 
as was previously mentioned, the calcium-argon baiid is 
weaker by about 70% than the argon-argon bond.B De- 
spite the relatively small size of our sample with respect 
to the number of existing structures, the heat capacities 
in Fig. H(a) exhibit a reasonably good agreement between 
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FIG. 8: Absorption spectra of the CaAr™ clusters corresponding to the heat capacities of Fig. ^. (a) n — 13; (b) n — 30; (c) 
n — 37; and (d) n = 54. The notations are the same as in Fig. 



simulation and the superposition method, the high peak 
near 55 K being the signature of the solidlike-liquidlike 
phase change. Melting is essentially caused by the mo- 
tion of argon atoms, and is relatively unsensitive to the 
location of calcium. This explains the shape and the evo- 
lution of the absorption spectrum in Fig. ^(b). The two 
peaks at 23720 and 24150 cm^^ become broader and van- 
ish between 25 K and 30 K, to be replaced by a major 
peak at 23600 cm~^. This transition is correlated with 
fluctuations in the location of the calcium atom over the 
cluster surface. On the other hand, melting at 55 K does 
not lead to any visible change in the absorption spectrum. 



D. Discussion 

The above results may not be immediately amenable 
to experimental comparison, because of the problem of 
size selection. Our simulation results show that com- 



plex isomerizations involving strong changes in the local 
environment of calcium could be, in principle, detected 
in experiments on size-selected clusters initially thermal- 
ized. However a typical sample of clusters in molecular 
beams is made of a broad size distribution, where the 
width and the center have values close to each other. The 
common feature of these clusters is the surface location 
of the calcium atom, and we may expect the absorption 
spectrum to be a possible probe of surface melting in 
these clusters, provided that the surface melting temper- 
ature does not significantly depend on cluster size. Even 
though calcium in bulk argoijuhas a very distinct, three- 
peak spectroscopic signature,llj it may be hard to detect 
volume melting in which calcium would be fully solvated. 
In fact, such van der Waals clusters evaporate atoms at 
temperature close to the melting point. From a statis- 
tical point of view, the weakly bonded calcium atom is 
likely to evaporate first. In addition the situations where 
calcium is fully solvated are less probable. Therefore the 



13 



4000 



3000 



o 

g_ 2000 

CO 

o 

■4— ' 

CO 
CD 

X 1000 



(a) Heat capacity 



Superposition 
' IVIonte Carlo 




20 



40 
T/K 



60 



80 



23000 



(b) Absorption spectra 




23500 



24000 



24500 



25000 



Wavelength (cm" ) 

FIG. 9: (a) Heat capacity and (b) absorption spectra of CaAri46- The notations are the same as in Fig. ^ 



true volume melting transition occuring with full solva- 
tion of calcium may hardly be observed. 

Provided that one could size-select these clusters, we 
have shown in the present work that spectroscopy can 
possibly probe their phase changes. More generally, in- 
tricate finite size effects in the thermodynamical behavior 
could be expected to be seen, as in CaAria, CaArao or 
CaAray. In these special clusters, coexisting isomers and 
multiple phase-likes induce complex variations of the ab- 
sorption properties with increasing temperature. 

Even with experimental setups involving selected sizes, 
another problem may come into play to hinder one from 
getting the complex spectroscopic information. All the 
above results are supposed to be at thermal equilibrium, 
without assuming how this equilibrium state had been 
reaching. In some cases, the relaxation time required 
for achieving equilibration may be very long because of 
multiple funnels in the energy landscape and kinetic trap- 
ping. This has been observed in the Arag cluster by Miller 
and coworkers,E3 where typical timescales of 10 s have 
been estimated for the interfunnel transitions at moder- 
ate temperatures. Relaxation to equilibrium and time- 
dependent absorption spectra will be the subject of the 
following article .Ej 



V. CONCLUSION 

In this paper, we have investigated some methodolog- 
ical and practical aspects of finite temperature spec- 
troscopy in chromophore-doped van dcr Waals clusters 
at thermal equilibrium. Using the harmonic superposi- 
tion approach improved with anharmonic perturbation 
corrections to the partition function, the basic thermo- 
dynamical observables were calculated from samples of 
local minima on the energy landscape. The photoabsorp- 
tion spectrum was obtained using the recent Gaussian 



theory by Wadi and PoUaknJ in the harmonic approxi- 
mation. These theories were tested and validated on se- 
lected CaAr„ clusters, in both the classical and quantum 
regimes, at low or moderate temperatures. Although an- 
harmonic corrections could also be incorporated to the 
various quantities (line shifts and widths) in this theory, 
we found them unnecessary in the present work because 
of the other underlying approximations. 

We have studied the possible spectroscopic signatures 
of the various isomerizations and phase changes in CaArg, 
CaArio, CaAris, CaArso, CaAr37, CaAr54, and CaAri46. 
The two smaller clusters show relatively simple isomer- 
izations leading to a similar spectroscopic signature with 
a clear three- to one-peak pattern. Both CaArao and 
CaAray display strong (non monotonic) finite-size effects 
and solid-solid transitions. In CaAria, melting involves 
the calcium atom taking place of one argon atom in the 
icosahedral shell. In contrast with most other sizes, it is 
associated with the appearance of a new absorption peak. 
In CaAr54 and CaAri4g, preliminary surface melting oc- 
curs at nearly 30 K, and is reflected in the absorption 
spectrum. In all these phenomena, the local environment 
of calcium changes rather abruptly, resulting usually in 
clear shifts in some absorption peaks. Volume melting 
towards a calcium atom fully solvated by argon has a po- 
tentially significant spectral signature, however calcium 
is less bound to argon than argon itself, hence it is likely 
to dissociate before volume melting occurs. One should 
mention that the main heat capacity peak indicates melt- 
ing of the argon cluster host, especially for large sizes. 
While the host cluster is liquidlike the calcium is still 
free, to some extent, to glide over the cluster surface. 
Being in this state, the rearrangements inside the argon 
cluster have little impact on the global spectroscopy of 
the chromophore, and the system can be expected to be- 
have as a pseudo diatomics Ca-Ar„|— This point will be 
examined further in the third paper E3 
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The techniques used here are vahd in both the quan- 
tum and classical regimes, which refers to the limit 
?i — > 0. They could be applied to other chromophore- 
doped van der Waals clusters, provided that potential 
energy surfaces are built first. In particular, neon clusters 
could be profitably investigated using the superposition 
method and the Gaussian theory of absorption. For he- 
lium, the harmonic approximation may break down, and 
even anharmonic corrections should be unsufficient. In 
this case, quantum Monte Carlo should be coupled with 
on-the-fly histogram accumulation to compute the spec- 
trum. More generally, the present methods are relevant 
to any statistical driven absorption spectrum in hetero- 
geneous or homogeneous polyatomic molecules. For in- 
stance, the thermal cffcicts-pn the absorption spectrum 
in small sodium clustery'lr'l could be investigated using 
appropriate modelsol£3 As concerns the approximations 
within the Gaussian theory, the most questionable point 
is the short-time expansion in the autocorrelation func- 
tion. This approximation could be_pvercome using the 
formal results of Yan and Mukamelo Unfortunately, this 
method involves a heavy numerical effort in the context of 
multiple isomers, and it still assumes harmonic forces. Fi- 
nally, the superposition approximation provides a frame- 
work to study time-dependent properties over long time 
scales not available to molecuiar dynamics, by solving 
appropriate master equations.LlI The application of this 
technique to photoabsorption spectra will be made in the 
next paper of this series E3 



phenomenological way 
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and we keep only terms containing coordinates up to the 
power 12, which contribute to corrections in the low- 
temperature classical partition function.Ea As in the clas- 
sical case, all odd products have a zero contribution, but 
it is no longer possible to sort the various nonzero terms 
according to their power in p. It appears more natu- 
ral to sort them according to the even power in coordi- 
nates. Following our previous notationsEj we write the 
anharmonic term 5V as a function of the normal mode 
coordinates Q as 
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The corrective terms can be expressed as averages on the 
harmonic ground state: 

Z{(3) w Zo(/3)[l + 74 + 76 + 78 + 710 + 712], (22) 

where 7^ labels the correction due to power i in coordi- 
nates. These terms are explicitely given by 
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APPENDIX: PERTURBATION EXPANSION OF 
THE PARTITION FUNCTION 



The ground state Hamiltonian is written in cartesian 
coordinates P,R as i?(P,R) = -ffa(P,R) + 5V{K), 
where Ha is the harmonic term and 5V the extra, anhar- 
monic contribution. The Taylor expansion of 5V near the 
minimum Ro contains products of at least 3 coordinates. 
The partition function Z{[3) = Tr exp^—PH) is now ex- 
panded in a semiclassical fashion, in order to yield at least 
the same terms as the corresponding classical expansion. 
This approximation is actually correct from the quantum 
point of view as long as we keep it to first order. How- 
ever, the rigorous perturbation treatment, which should 
be performed with path integrals, will not allow us to 
obtain the simple corrections that can be compared with 
the classical expressions of Ref. pH. Thus we write in a 
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The above expressions can be considerably simplified by 
factorizing the average products {Qi x • • • x Qj). This is 
made possible by Wick's theorem, which holds for general 
bosonic ground states: 
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The two diagrams giving 75 are 
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It is then more convenient to write the 7i's as diagrams, 

(k) 

in which each A:— vertex represents one of the V\_^ ...^^ term. 
For instance we have 



74 



with 



00 = ErS(Q?>(Q?) 



(30) 



(31) 



76 




^ *^+^ (32) 



The other perturbative terms 7g, 710, and 712 contain a 
larger number of diagrams, including non connected ones. 
The fewest diagrams are obtained in In Z, which still has 
46 different contributions corresponding to Eq. (22): 
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It is important to notice that this perturbative expan- 
sion includes all-IB terms up to second-order in the clas- 
sical expansion,t£l and that 28 remaining terms would 
contribute classically to 3rd and higher orders in T". 
These terms correspond to the T^^) x T^^), T'^^I x T^^), 
r(3) X r(4) X r(5), and r^^) x r^") x r^^) products, but 
cannot be neglected a priori with respect to the other 



(33) 

I 

corrections in the quantum regime. 

In the superposition method framework, these pertur- 
bative corrections act on the individual partition func- 
tions, where they are expected to give a more accurate 
estimate of the relative weights of the isomers as a func- 
tion of temperature, as well as more accurate thermody- 
namical observables. 
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